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Abstract 

The nonlinear viscoplastic behavior of fibrous periodic composites is analyzed by discretiz- 
ing the unit cell into triangular subvolumes. A set of these subvolumes can be configured by the 
analyst to construct a representation for the unit cell of a periodic 

the loading history the total strain increment at any point is governed by an integral equat 

entire composite. A Fourier series approximation allows the incremental 
^esse^d^ stains to be deternLd within a unit cell of the periodic lattice. The nonlinearity 
arising from the viscoplastic behavior of the constituent materials comprising the composite 
treated as a fictitious body force in the governing integral equation. Specific numerical exam- 
ples showing the stress distributions in the unit cell of a fibrous tungsten/copper metal matrix 
composite iLder viscoplastic loading conditions are given. The stress distribution result g 
the unit cell when the composite material is subjected to an overall transverse stress loading 
history perpendicular to the fibers is found to be highly heterogeneous, and typical homo 0 
nization techniques based on treating the stress and strain distributions witlun the constituent 
phases as homogeneous result in large errors under inelastic loading conditio • 


1 Introduction 

Structural analysis of component hardware fabricated from composite materials requires the ma- 
terial to be considered as homogeneous. Methods for developing homogenized constitutive models 
for the composite rely on volume averaging the local stress-strain response in a small representa- 
tive volume element of the composite material. Most methods approximate the local stress and 
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strain fields in the representative volume element as being homogeneous in each phase. If higher 
accuracy is required finite element representations of the local response can be obtained, but such 
finite element analyses would have to be carried out at each integration point in the global finite 
element analysis of the component at each step of the loading history. Such an approach leads to 
prohibitively long computation times. 

In order to increase the accuracy of the local stress and strain field representations, and the 
efficiency in obtaining them, methods have been developed which can be embedded in the finite 
element calculations at each integration point. In previous papers 1 ’ 2 we developed a procedure for 
analyzing the nonlinear deformation behavior of periodic composites by means of Fourier s series 
and Green’s functions. These methods are built upon work carried out in references 4—9 and 10-18, 
respectively. Both of these complementary approaches involve the solution of an integral equation 
for the total strain increment, in which the nonlinearity arising from viscoplasticity is treated as a 
fictitious body force. The loading history imposed on the composite is divided into discrete load 
steps, and the integral equation for the displacement increment in each load step is derived from 
Navier’s equation when subject to incremental viscoplastic body forces. In the references cited 1 ’ 2 , 
both approaches were shown to be equivalent to one another, and the Green’s function formulation 
can be derived from the Fourier series representation by means of a Poisson sum technique. In 
another paper 3 the Fourier series approach was used to analyze the elastic behavior of a tungsten 
fiber/copper matrix (W/Cu) composite which is being considered as a substitute liner material 
to increase the strength and improve the durability of the combustion chambers in cryogenic ally- 
fuelled liquid propellant rocket engines. This analysis was carried out by discretizing the unit cell 
into rectangular subvolumes. In this paper we introduce a triangular discretization of the unit cell 
and perform a nonlinear viscoplastic analysis of periodic (W/Cu) composites. The rectangular and 
triangular subvolumes are the first two element geometries of a library of geometries, similar to a 
library of finite element geometries, that a user might wish to employ when performing a composite 
analysis. 

In the following sections we develop numerical techniques which can be used to analyze the 
local stress and strain fields within the unit cell of a periodic composite. In an actual finite element 
calculation the unit cell would be discretized into a small number of subvolumes, and the local strain 
field values would be homogenized to provide the finite element code with the overall homogenized 
constitutive response at each integration point in the structure. Post-processing of the finite element 
results can then be carried out using a much finer discretization of the unit cell with the proposed 
numerical techniques. 
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2 Constitutive Equations 

2.1 Hooke’s Law 

The form of Hooke’s law suitable for infinitesimal strains can be written in the incremental form 

Ac,, (r) = Dijki (r) (Ae£ ( r ) _ Ac*, (r)) , (!) 

in which Ach (r) is the strain increment representing the deviation from isothermal elastic behavior, 

A CM (r) = Ae\, (r) + a*, (r) AT (r) , ( 2 ) 

where at the point in the unit cell with position vector r, D,, tl (r) is the elasticity tensor, Ael, (r) 
is the inelastic strain increment, and au (r) AT (r) is the thermal strain increment. The inelastic 
strain increment can be computed explicitly at the point r since the stress is assumed to be known as 
a function of position r at the beginning of the increment. The only unknown quantity in equatron 

(1) is then the total strain increment, A e& (r). 

2.2 Evaluation of Total Strain Increment 

The elasticity tensor at any point r in the composite material may be written in the form 

d«hW = BJjH + MWO, (3; 


where 

SDiju (r) = ’) (r) (d\,u - D T,ti) ■ 

In this relationship t?(r) = 1 in the fiber and 0 (r) = 0 in the matrix, with D\ jU denoting the 
elasticity tensor of the fiber and Dg u that of the matrix. We then fad 1 ' 3 that in the Fourier 
series approach the total strain increment in the unit cell of a three dimensional periodic lattice is 

determined by solving the integral equation, 

1 ±cx> i 

Ae£j ( r ) = ^ £ kl + 77 X) X X 9klmn (C) X 

7l p = 0 

x JJJ {D™ nrs Ac rs (r') - 6D mnT3 (r') [Ae* (r') - Ac rj (r')] } dV{ r') , (5) 
where the fourth rank tensor g k i mn (C) is § iven b y 


gklmn (0 = 2 (CnClM mk (C) + CnCfc^m/ (O) > 
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in which the Christoffel stiffness tensor M tJ (£), with inverse M {j 1 «), is defined 11 ’ 12 by the relation 


M tj «) = D%'C, 


(7) 


with C P = = £ P /£ being a unit vector in the direction of the Fourier wave vector £, and 

4 = denoting the magnitude of the vector £. In equation (5) the sum is taken over integer 

values in which 


6 = 


27rni 


6 = 


2tTU2 


& = 


27T713 


(8) 


L J ^ r ) t 

1 xv2 

and L 1? L 2 , L 3 are the dimensions of the unit periodic cell in the x u x 2 , £3 directions, so that the 
volume of the unit cell is given by the relation, V c = LiL^L^. The values of ni, 7i 2 , 7i 3 are given by 


n p = 0, ±1, ±2, ±3, . . . , etc., for p= 1,2,3 (9) 

and the prime on the triple summation signs indicates that the term with ni = n 2 = n 3 = 0 is 
excluded from the sum. 


2.3 Evaluation of Viscoplastic Strain Increment 

In solving the preceding integral equation for the total strain increment, AeJ (r), it is first necessary 
to determine the inelastic strain increment, Aeh (r), and thence the deviation strain increment, 
A c k i (r), in equation (2). 

For this purpose we use a simple viscoplastic constitutive relationship whose governing equations 
at a temperature of T° K are given in the form 19 : 


PI 

6(T) 

l|S|| 

S{j 



r 

exp 

1 

50 

1. 1 


> 

for T t <T <T m 

\ 

exp 

q 

RT 

HIM 

for 0 < T < T t 

J j 



&ij 3 


( 10 ) 

( 11 ) 


( 12 ) 

(13) 

(14) 


Material constants for this viscoplastic constitutive formulation 19 for copper have the values: 
A = 2 x 10 7 per second; C = 13 MPa; n = 4.5; Q = 2xl0 5 Joules/mol.; R = 8.314 Joules/ (mol. 0 K); 
T m = 1356° K; and T t = \T m . 
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The preceding equations represent a simple Norton creep law where no account is taken o e 
history dependence of the inelastic deformation behavior. This is adequate for the type of loadrng 
considered in this paper, but for more complicated loading conditions it is necessary to use some 
form of unified viscoplastic constitutive formulation'’ where history dependence rs accounted 

with internal state variables. 

3 Numerical Solution of the Nonlinear Integral Equation 

3.1 Solution Under Overall Strain Control 

For two-dimensional fibrous composites the integral equation for &*(») reduces to the form 

A.S(r) = M. + jEEWC )JJe^{D^A^)- 

- 6D m „, (r') [A4 (r'J - Ac„ (r 1 )]} dS(r') , ( 15 ) 

where A c = L.Ls is the area of the unit cell. Nemat-Nasser and his colleagues*-’ and the authors'- 3 
have demonstrated that good accuracy can be achieved by dividing the unit cell into a number 
of subvolumes, and by approximating the strain increment A* M in the /3th subvolume integral 

with its average value in the subvolume, viz., 


AeJ, (r') = AeJ/ = j~ jj dS ^ ’ 


( 16 ) 


where A s is the cross-sectional area of the subvolume. 

Let there be N subvolumes in the unit cell, with M subvolumes in the fiber and N - M 
subvolumes in the matrix. Then the preceding integral equation, with a piecewise constant strain 
approximation in the subvolume integral, can be written as 

A.S(r) = A sU±ti^S«™W ( 'JJ e - <(T ' iS{r ') X 

"C /3=1 «p=0 A s 

x {D”„ rs A£ - [AeJ/ - Ac?,]} , ( 17 ) 

where = D’ mnr , - D: or 0, according as the subvolume 0 is in the fiber or matrix, 

respectively. 

If we use Nemat-Nasser’s notation and write 


1 


and denote 


( 19 ) 


fa _ 

1 A e 


as the volume fraction of the ath subvolume, then the preceding equation may be written as 


N ±oo 


IV _L. r _ -i n 

A4 (r) = Ae° t , + EEE' 9k h» «) e*'fQ e (- «) [Ae™ - A<£] } . (20) 

0=1 n p =0 


We may now volume average equation (20) over the ath subvolume to obtain 


N 


M 


A 4“ = A4 + E - E 

7=1 0=1 


where 


and 


±00 


SS?.= EE (0 sr&„,<y(t) </(-«) 


±co 


<„ = s£ + E E' 9 h~ (0 </(-«)• 


( 21 ) 


( 22 ) 


(23) 


n p =0 


The matrix tensors Alf r3 and S$ rs are akin to Eshelby’s 10 - 11 tensor for an ellipsoidal inclusion, but 
account for the interaction of the fiber with its neighbors in the infinite lattice. It may be noted that 
the last term in equation (21) is summed only over the subvolumes in the fiber, where (1 < (3 < M), 
since 8D&. = 0 if the /5th subvolume resides in the matrix, so that 


S?P = 0 for M < (3 <N. 


(24) 


Thus only M unknowns (associated with the subvolumes in the fiber) are involved in equation (21). 
When this relation is assembled columnwise for each subvolume a, the solution to the 6M x 6M 
system of equations can be obtained by L U decomposition. We then have to solve the system of 
equations 

M 

E B ^A^ = m for a = 1.2 (25) 


0=1 


where 


Bg, = + fSg,, (26) 

with Ikirs = (<Sfcr<5ii + SkaSir )/ 2 denoting the fourth rank identity tensor, <5 a/3 the matrix Kronecker 
delta, and 

A b° kl = A 4 + E PAHAcl,. (27) 

7=1 
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Equation (21) can then be nsed with the known values of A «JT in the fiber snbvolnmes, where 
(1 < p < M), to compute the values of A eft in the matrix subvolumes, where (M < a < W). 

An explicit solution for the total strain increment can be found from equation (25) if the total 
strain increment is small Under these circumstances the inelastic strain increment can be computed 
explicitly from equation (10) in the Euler forward difference form, 


m 


IiL At. 
IISII 


(28) 


In this form the inelastic strain increment is independent of the total strain increment and the right 
hand side vector Ai>& in equation (25) can be evaluated from equations (2), (27) and (28). However, 
if an implicit scheme such as backward Euler integration (or a forward difference submcrementation 
method) is used to evaluate the inelastic strain increment, then As?, will depend on As?, and 
equation (25) must be used iteratively to solve for A eh , ie - 


E^{ A ^L. = Ai, “({ Ae «h) 


for 


( a = 1,2,... ,M 1 
\ 7 = 1,2,. ..,1V | 


(29) 


where A is the iteration number* 

Once the values of AeJ° have been determined, further resolution can be obtained, if required, 
from equation (20). By substituting different r values in this equation we may evaluate AeJ(r) 
at any point within a subvolume. However, the Fourier series solution in (20) will exhibit Gibbs’ 
phenomenon, in which the solution oscillates spatially with r near the fiber/matrix interface. These 
oscillations have a period 20 of L,/Af, for 7 = 1,2 where L, is the length of the unit cell and JV, is 
the number of terms in the Fourier summation corresponding to the 7 th direction. It is therefore 
appropriate to average equation (20) over a small rectangular area about the point r. This results 
in the term e‘« r being replaced with the Laue interference integral, <?“({), for a small rectangular 
subvolume 3 of dimensions Li/Ni x L 2 /N 2 - 


3.2 Solution Under Overall Stress Control 

In the preceding section the equations governing the total strain increment distribution in the unit 
cell were derived under the assumption that the composite was subjected to a uniform overall 
strain increment given by A e%. When a uniform overall stress increment, AoJ, is applied to the 

composite, it is necessary to obtain As 3 - as a function of A <7°-. 

The overall stress and strain increments are equated to the volume averages of A<r* (r) and 
Ae T( r ) over the unit periodic cell. Hence, for fibrous composites, equation (1) can be volume 
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averaged over the unit cell to give 

Aa ?■ = ^JJ Dau (r) (Ae£ (r) - Ac« (r>) d5(r). 


( 30 ) 


When equation (3) is substituted into the preceding equation and we replace A s kl (r) and A c k i (r) 
with their discretized subvolume counterparts, we obtain 

(31) 


in which 


and 


Aa°. = A-^Aefc + + A7 *’ 


Aft, = £ /* [(A" - A”) Aejf + 2 (/ - /i") tof\ 


(32) 


0=1 


AD; -Ef (A^AcJ, + 2^Acj) . (33) 

In the preceding equations isotropic relations have been assumed for D? jkl and 6Df jU in the form 

D iikl = ^ (hkSjl + tiirfjk) ( 34 ) 


and 


6D? jkl = (\ 0 - A m ) SijSki + (/ - V m ) (M;/ + SaSjk ) , ( 35 ) 

where the superscript 0 = f or m, according as the /3th subvolume is in the fiber or matrix, respec- 


tively. 

Equation (31) can be inverted to give 


„ Ag?, - Aft, - ad, A^ 

A£ ij - 2 n” ” 2fi m (3A m + 2/i“) 


(A<4 - Ai?* fc - AT fcfc ) • (36) 


This may now be substituted into the integral equation (21), and when all the terms containing 
A £ Ta are taken over to the left hand side, the system of equations takes the form 


M 


where 


£ C&At™ = Aaa for a = 1,2 


ctL = « a ‘ s 4ir, + /' , sat + ^:{(A 9 -A")M„ + 2(/-^ m )4 I r.}- 


(37) 


f 


2/x m (3A m + 2/i m ) 


{3(A £ -A”)+2(/-m”)}W, 


(38) 
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and 


(39) 


Aa£ ; — 


A a®/ - A Tki 


a _ 


2/i 


m 


{„ £ . (A<r° - ATpp) + E r^S.A<?.- 

° kl 2/z m (3A m + 2^ m ) V m J 7 =i 


When equation (37) is solved for A e*f for all the subvolumes in the fiber, Aft,- may be evaluated 
from (32) and substituted into (36) to obtain the overall strain increment, AeJ. Equation (21) is 
then available to determine Ae£ a in every subvolume in the unit ceU. 

3.3 Matrix Assembly of Discretized Integral Equation 

The solution of the integral equation in (25) requires the evaluation of the matrix tensors S& and 
AtL from equations (22) and (23). For isotropic constituents the matrix tensor Sg rt involves the 
evaluation of g klmn «) ^i nT3 , which may be written from (6) and (7) in the form 

f - \™\ ( A m 


Qklmn (C) ^mnrj 


A E 


X m + 2fj, m i 


SrsCkCl + 


+ ^ {SrkCsO + fas CrO + hrCaCk + ^IsCrCk) 

V A Am + 2 ^ / 

where X 13 , // and A m , /j m are the Lame constants for the /3th subvolume and the matrix, and 
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Ci = 


for z = 1,2- 


(40) 


(41) 


y(27mi/Li) 2 + 

with 6 = ixni/Lu (2 = 2tra 2 /L 2 . The assembly of the matrix tensor requires the evaluation 
of 9Hm , «) IC.„. and this is obtained by replacing A* - A" and / - in equation (40) with A” 

and /i m , respectively. 

In the preceding derivation the isotropic form of equation (7), 


with inverse 


Mij (C) = /*“%• + ( A “ + M 1 ") CtCji 


c Am i 


(42) 


(43) 


obtained by substituting (34) into (7), is used. Substitution of (43) into (6) together with the use 

of (35) then yields the relationship given in (40). 

The remaining factor required for assembling the matrix tensors S&, and A&. is the Laue 
interference product, Q*({) </(-«), with «“(«) defined by equation (18) in which A a denotes the 

area of the ath triangular subvolume. 
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Let the vertices of the ath triangular subvolume have coordinates (x?,yf), (xg,!/?)’ ( x 3iVz)- 
We map the Cartesian coordinate system ( x,y ) into a new system (r,s), in which the vertices of 
the mapped triangle are at the points (0,0), (1,0), (0,1). This mapping is accomplished with the 

relationships 

x = x? + (xf -x?)r + (xg -x%)s ( 44 ) 


and 


y = Vi + (y? - Vi) T + (2/3 - Vi ) s • ( 45 ^ 

The area transforms according to the usual Jacobian relation dxdy = J (r, s ) dr ds, in which 


9(x,y) _ 

dx dy 
dr dr 


x? - x? y£ - y? 

5(r,s) 

dx dy 


xg “ *1 2/3 “ 2/? 

ds ds 




(46) 


(47) 


J ( r , s) = 

= (xg - x?)(yf - y t) - ( x t - x t)(y 2 - vi) = 

and where the last equality is obtained by means of the identity, 

JS d xdy = £ o /£' J (r, >) ir is = = A.. 

A a 

The Laue interference integral can now be written as 

Q“(£) = J_ /"/" e i( ^ xHiy) dx dy 

A.Q 

_ 1 f 1 f 1-3 p iUi[xf+(x9;-xf)r+(xf-xf) 3 }+^[yf+{y?-y?y+{yf-y?)‘>]} dr ds 

- T a LoJr=0 5 ( r - S ) 

2 f 1 / 1 "% 1 {« 1 K + ( I ?- I “) r+ ( I r-*f>]+€a[»f+(»?-»fW»»‘-J'f)']}dr(i5. (48) 

J3 = 0 Jt = 0 

After the integration is performed, the relations for Q Q (£) and can be separated into their 

real and imaginary components, yielding 

$“(*) = + (49) 


and 

Q 0 (-Z) = Q?(£ i,6) - *QS(4i»e*)- ( 5 °) 

In evaluating the Laue product Q a ($)^(“0» we need retain only the real part, since the imaginary 
part vanishes on summing over m,n 2 in the Fourier expansion. Accordingly, we may write 

Q a (0<?*(-0 = Qtitu&Qi&’b) + 6)«?(ei. 6)- ( 51 ) 
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On integration, equation (48) yields the relations, 

cos (frxg + £ 23 / 2 ) 


Q?(€ 1.6) = 2 


and 


Q2(6,6) = 2 


(^l( x l — x !z) + 6(l/f — 2/2 )) (6( x 2 x ?) + &(y2 2/3)) 

COS (|jXj* + & y? ) 

( 6 ( x ? - x ?) + 6(2/1 - 2/2)) ( 6 ( x “ - x ?) + 6 (y? _ 2/3)) 

cos (gixg + 62/3) 

( 6 ( x ? - x ?) + 6 (yf - 2/3“)) ( 6 ( x ? - x ?) + 6 (y? - 2/3)) 


sin (fraff + 6^2 ) 

(6( x “ - x 2 ) + 6(y? - 2/2)) (6( x ? - x 3) + 6(2/2 - VS )) 
sin (6 x i + 62/1) 

(6( x “ - x ?) + 6(y“ - y“)) (6( x “ - x s) + 6(y? - y? )) 

sin (|ix| + 6y“) 

(6( x “ - x ?) + 6(yf - y?)) (6( x ? - x s) + 6(y? - y“)) . 


( 52 ) 


(53) 


The relations for Q? (6, 6) and Q\ (6, 6) are the same as the preceding relations with 0 substituted 
for a. 

In the case of rectangular subvolumes we have shown 3 that the Laue interference integral con- 
tains products of the function sinx/x and it is necessary to attend to the limiting case where the 
denominator x -> 0 and sinx/x - 1 as x - 0. For triangular subvolumes the denominators in 
the Laue interference integrals in (52) and (53) vanish under nine separate circumstances, and the 
limiting form of the integrals, Q?( 6,6) and Q?( 6,6), are given in the Appendix. The confluent 
case where & and & are both zero is of no concern, since this case is- excluded from the Fourier 

summation. 

In the limiting forms given in the Appendix, it is assumed that the denominators are exactly zero. 
When, for example, the denominators are nonzero but small, it is necessary to expand equations 
(52) and (53) into a Taylor series in ascending powers of the small quantities for the nine cases. 
This can easily be achieved with the aid of a symbolic algebra software package. 

4 Numerical Examples 

Bahei-El-Din, Dvorak and their colleagues 21 ’ 22 have given cogent arguments which show that it is 
inappropriate to assume that the stress and strain increments in the unit cell are spatially uniform in 
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each phase of the composite. They found that when the true (within discretization error) stress and 
strain increments in the unit cell were calculated throughout the loading history and the results were 
volume averaged to produce the overall macroscopic response, the results were significantly different 
from those (e.g. self-consistent and Mori-Tanaka schemes) in which the overall macroscopic response 
was found by assuming spatially constant fields in each phase. This was found to be especially true 
in composites loaded into the plastic region where the constitutive material properties vary at each 
point in the unit cell in accordance with the spatially varying field histories. In the preceding studies 
the true field histories in the unit cell were obtained with the use of a periodic hexagonal array 
model 23 in which periodic boundary conditions were applied to a finite element discretization of the 

unit cell. 

Aboudi and his colleagues 24 " 26 have also examined the field histories within the unit cell of 
a periodic composite by discretizing the unit cell into rectangular subvolumes and by enforcing 
periodicity through continuity of displacements and surface tractions at the unit cell boundary. 

Both of the preceding methods employ a field theory approach in which the deformation behavior 
within a subvolume depends only on the displacements and tractions at its surface. Here, we study 
the field histories by discretizing the unit cell into triangular subvolumes, and periodicity is enforced 
naturally by expanding the field histories into Fourier series. Since an integral equation scheme is 
adopted, the method is an “action at a distance” approach in which the deformation behavior within 
a subvolume depends on the field histories in all the other subvolumes comprising the composite 
material. The Fourier series method essentially restricts the interaction between the subvolumes 
to those within the unit cell, whilst the Green’s function approach 1 " 3 ’ 28 requires the interaction to 
be evaluated between all subvolumes in the composite, though the interaction is strong only for a 
subvolume which interacts with the subvolumes in its own and neighboring unit cells in the penodic 

lattice. 

The size of the matrix which must be solved by LU decomposition for a fiber with M subvolumes 
is 6 M x 6M. If the fiber has a square planform with n subvolumes along each side, then the size of 
the matrix increases with the number of subvolumes along a side according to 36n 4 . For a three- 
dimensional cuboidal inclusion the size of the matrix increases as 36n 6 , so that a three-dimensional 
problem has a matrix size n 2 times as big as that for a corresponding two-dimensional problem. 

In order to demonstrate that the triangular subvolume approach gives reliable results, we first 
calculate the result of elastically loading a fibrous metal matrix composite. Figure 1 shows the 
transverse stresses, a n , in each subvolume within the unit penodic cell, when a W/Cu fibrous 
composite is loaded in the transverse direction with an overall uniform stress of ah = 1000 kPa. 
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The tungsten fiber is assnmed to be of sqnare planform and occapies a volume fraction of / 
9/49 = 0 184 in the unit cell of the composite, with Ew = 395 GPa, fw - 0.28, £c« 

GPa and , c „ = 0.34. Fig. 1(a) presents a numerical tabulation of the constant-valued stresses 

for each of the 49 square subvolumes, using the Laue interference product, <3 ‘(i)Q (-«• tor 
rectangular subvolumes presented in a previous paper*. The corresponding results when the unit 
cell is discretized with triangular subvolumes is presented in Fig. 1(b), where it may be seen 
the ceU is discretized with the hypotenuse of each triangle biased in the northeast direction. is 
causes slight asymmetries in the resultant elastic stress distribution within the unit cell. However, 
when the results from the two triangular subvolnmes in Fig. 1(b)— into which each rectangular 
subvolume in Fig. 1(a) is broken-are averaged, the results in the two figures agree very close y. 
These results were obtained by summing the Fourier series to n, = ±200 for p - 1,2 in equations 
(22) and (23). The results do not differ significantly from those obtained by summing the series to 


P We now consider the W/Cu composite as an assemblage of unit cells in which the tungsten fiber 
is assnmed to be elastic whilst the copper matrix is assumed to be a ductile viscoplastic material 
wh»e inelastic response is governed by equations (10)-(14). The triangular discretization of the 
unit cell is shown in Fig. 1(b). Another discretization, in which the hypotenuse of each triangle 
is biased in the northwest direction, was also used and the subvolume stresses were averaged to 

eliminate asymmetries due to meshing. 

The composite is loaded in the transverse direction under overall stress control at a stress rate of 
d« = 10 MPa per second to a maximum overall stress of = 100 MPa at a temperature of 814“ K 
0 ‘ 0 .6T„. It is then allowed to creep at this overall stress for 3600 seconds, and is then unloaded to 
an overall stress of <ft = 0 MPa at an overall stress rate of oft = -10 MPa per second. All other 


components of other than oft are taken to be zero. 

Figure 2(a) shows the transverse stresses, on, within the unit cell at the end of the mitral 
second loading when the overall stress is oft = 100 MPa. Stresses are rounded to the nearest 
integer. The numerical tabulation in Fig. 2(a) is presented as a topological map in Fig. 2(b), where 
the height represents the magnitude of the transverse stress concentration factors. The transverse 
stress distribution forms a ridged valley where the stress gradient parallel to the loading direction is 
relatively small, and where the stress is higher in the ridges near the north and south fiber/matrix 
interfaces and is smaller in the valley which runs parallel to the loading direction in the center of 
the fiber. The formation of a ridged valley is due to the assumed sqnare shape of the tungsten fiber, 
and to a lesser extent, to the interaction between the neighboring fibers in the periodic lattice. A 
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similar ridged valley is observed when a cuboidal fiber embedded in an infinite matrix experiences 
a uniform eigenstrain 11 . 

' Figure 3(a) shows a tabulation of the local hydrostatic stress field, (an + a 22 + cr 33 )/3, for the 
same transverse loading condition. High hydrostatic stresses occur at the fiber/matrix interface 
perpendicular to the direction of the loading axis, and Fig. 3(b) shows the results as a topological 

map. 

The transverse stresses within the unit cell after one hour of creep are shown in Figs 4(a) and 
4(b). High stress gradients are present at the fiber/matrix interface and creep has the effect of 
increasing the heterogeneity within the unit cell in both the fiber and the matrix phases. On 
examining the hydrostatic stresses in Figs 5(a) and 5(b), we observe that some parts of the fiber 
are in hydrostatic tension whilst others are in hydrostatic compression. The matrix is in a state of 
hydrostatic tension with very large values observed at the fiber/matrix interface perpendicular to 
the loading axis. 

After the composite is unloaded to a zero overall stress state, we see from Figs 6(a) and 6(b) that a 
high degree of residual transverse stress heterogeneity exists within the unit cell with relatively large 
tensile stresses in the ridges and compressive stresses in the central ridged valley of the fiber. The 
hydrostatic stresses in Figs 7(a) and 7(b) show a similar heterogeneous trend, with large compressive 
hydrostatic stresses in the fiber’s central valley and smaller compressive hydrostatic stresses in the 
ridges. In the copper matrix, the tensile hydrostatic stress adjacent to the fiber/matrix interface 
perpendicular to the loading axis is about as large as that which would obtain by tranversely loading 
a unif orm matrix material to 100 MPa. 

As Dvorak has aptly remarked 22 , these results serve as a salutory reminder that large errors 
can be. obtained in homogenization techniques which employ the assumption that stress and strain 
distributions are homogeneous in each constituent phase of the composite under inelastic loading 
conditions. 

5 Conclusions 

This paper presents one of three known methods for analyzing the heterogeneous deformation 
behavior in the unit cell of an infinite periodic lattice when subjected to a uniform applied stress 
or strain loading increment. An integral equation approach is used to determine the total strain 
increment at any point in the unit cell. The integral equation is solved by a Fourier series expansion 
of the field variables and their spatial variations within the unit cell are approximated by a piecewise 
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constant distribution. Other approaches which have been considered in the literature are A ou i s 
method of cells 21 - 26 , the Suite element method with periodic boundary conditions 2 '- 22 , and boundary 
element methods 22 . The inelastic stress distribution in the unit cell of a W/Cu metal matrix 
composite is evaluated when subjected to an overall transverse loading history imposed at a constant 
overall stress rate. The stress distribution in the unit cell is found to be very heterogeneous with 
the heterogeneity increasing with continued inelastic deformation in both the fully loaded and the 

unloaded condition. 
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7 Appendix 


Case 1 : y| = y? and 6=0 

2 (cos (6y?) - cos (6y? ) + 6 (y“ - y?) sin (6y? )) 
Q?(6>6) = - 6 2 (y? - y?) 2 

and 2 (sin (&y?) - sin (6y?) ~ 6 (y? - y? ) cos (&y?)) 

Q?(6,C2) = CKyf -y?) 2 


Case 2: y“ = y? an( ^ Ci = 0 

2 (cos (&y?) - cos (6yf) - 6 (y? - y? ) sin (6y? )) 
= 6 2 (y? - y“) 2 

^ ^ ^ 2 (sin (6 y? ) _ sin (6y? ) + 6 (y“ - ySpcos^yff)) 

Q2(€i*4a) - ' Wjj? - y?) 2 


Case 3 : y? = y? and 6 = 0 

2(cos (£ 2 y?) - cos (6y?) + & (y? - y?) sin (6y? )) 

and 2 (sin (&l/f) ~ sin (&y? ) — & (yf — Vs) cos (bvf )) 

= ■ <?(»? - ss) 2 


Case 4 : x? = x? and 6 = 0 

Q? (&»&) = 

and 

Q?(6,-&) = 


2 (cos (6a?) - cos (6a?) + 6 (x? - x?) sin (fra?)) 
£i(a? - *^2 ) 2 

2 (sin (^a?) - sin (6a?) - 6 (a? - a?) cos (6a?)) 
6(a? - x?) 2 


Case 5 : x? = x? and 6 = 0 


2 (cos(?iX?) - cos($iif) - 6 (if - 15) sin (fcif )) 

Qitiub) = g(*f - 15) 2 


( 54 ) 

( 55 ) 

( 56 ) 

( 57 ) 

( 58 ) 

( 59 ) 

( 60 ) 
( 61 ) 

( 62 ) 
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and 


2 (sin (frag) - sin (fra?) + ft (x? - sg)cos(ftsg)) 

Q?(fti6) = - x?) 2 


( 63 ) 
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Case 6: x“ = x“ and £2 — 0 


<W(£i,6) = 

and 

<?f(£x, 6 ) = 

Case 7 : Ci (gf - *2 ) + 6 (yf - 

Qf(Cx> 6 ) = 2 


and 

Qf(Ci,&) = 


Case 8: Ci (xf _ gf) + £2 ( 2/2 


2 (cos (Cigf) - cos (Cigf) + (gf ~ xg) sin(Cixf)) 

C?(gf - gf i 2 

2 (sin (Cigf) - sin (Cigf) - Ci ( xf - xf ) cos (Cig?)) 
(i< x i - x °) 2 

• vs) = 0 


( 64 ) 


( 65 ) 


Ci(gf (2/? - 2/3 ) - *2 (yf - yf) - *3 (»? - ))' 


(zf - xf) 2 cos 


| &(gf (yf-yf)-yf(*f -gf)) - 


(g? - gf) 


Ci (gf (yf - yf ) - *2 (yf - yf ) - *3 (yf - yf ))' 


+ 


(xf - xf ) sin 


f 


, C 2 ( g?y2 - x fyf) \ 
(gf-gf) ) 


6(*f (t/S - vs) - z °- (VS - VS) - x °(vS- VS)) 


( 66 ) 


Ci(gf (yf - yf ) - X Z (yf ~ tf) - x “ (yf - yf)) - 
C2 (gf (yf-y^-yf^f-^f))' 


(xf - x£) 2 sin 


(xf gf ) 


(l (if (vs - vs) - x % (vS - vt) - X S (vS-vS)) 

' . . . . {b(*SVS - *SVS)\ N 

(if - if) cos ( (xf _ x j— ) 

&(*? (vs-vs)- x s (vs - vs) - x %(vs- vs)) 


( 67 ) 


-VS) = 0 
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«f(6,6) = 2 


/ a a \2 ( - *3#) ^ 

(x 2 x 3 ) COS y ^ gg) J 


il (sf (t/f - y?) - *2 (yf - yf ) - *f (y? - tf ))' 
'6 (if (y 2 - y? ) - y? - ^3 )) ' 


(xf - if) 2 cos 


(*2 “ X f ) 


g(*f to? - vt) - x t (vt -vt )- x J to? - »>))' 


+ 



(if - if 

' ) sin ^ 

&(ifyf -ifyf)) 
(if - if) J 


6 1 

(i°( 

vt - vt) 

l-s?l 

(yf - yf ) - *f 1 

(vt 

: - yf )j 


+ 


(68) 


and 


/ 


QfUi,6) = 2 


(s 


-if) 2 sin ^ 


6(^2 2/3 ~ S?l/f A 

(l? - g?) / 


& 2 (if (2/2 - yf) - *s (yf - 2/?) - *? ( y % - ti))‘ 


+ 


(x% - £?) sin 


/6(g? (2/2° - 2/3) -yf (*2-33))’ 


(^2 ^3) 


& 2 (i? (2/2* - yf) - if (yf - yf) - ®3 (yf - yf ))' 

, a tt N / £2(ifyf ~ x 3 yf ) \ \ 

(x 2 X 3 )C0S^ (ga-ga) J 

& (xf (yf - y 3 a ) - if (yf - y? ) - if (yf - 2/?)) 


(69) 


Case 9: & (if - if) + 6 (yf - y? ) = 0 


Qf(Ci,6) = 2 


' (if - 1 ? 

) 2 cos I 

( 6 (ifyf - ifyf)\ 
l (if - sf) ) 

& 2 (zf (yf-yf) 

(xf - if ) 2 cos j 

-if 1 
1 £ 2 (if 

( 2 /f-yf)-sf (yf-yf))' 
(yf - yf ) - yf (if - *?)) 

v (If - 1?) 


$ (if (yf - y?) - if (yf - y? ) - if (yf - yf ))' 


+ 
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(a) Unit cell of 49 squares subvolumes 
with 9 squares ub volumes In the center 
embedded in 40 square subvolumes. 



(b) Unit cell of 98 trianglar subvolumes with 18 tri- 
angular subvolumes in the center embedded In 80 
square subvolumes. 


Figure 1 .—Transverse stress concentration, a-j-j , for an applied 
stress oT cr^ = 1000 kPa at an overall stressed rate of = 
1 000 MPa per second. Each unit cell is embedded In a doubly- 
periodic array of identical cells. The subvolumes in the center 
of the unit cell represent a tungsten fiber of square planform 
embedded in subvolumes representing a copper matrix. 
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(a) Transverse stress concentration. 




Figure 2.—' Transverse stress concentration, 1 , for an applied 
stress of or^* 100 MPa at an overall loading rate of =* 

1 0 MPa per second, and topological plot of the numerical 
stresses. 
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(a) Hydrostatic stress concentration. 



(b) Topological plot. 


Figure 3. — Hydrostatic stress concentration, ( a 11 + a 2 2 + <*33)^. 
for an applied stress of 0^= 100 MPa per second applied at an 
overall loading rate of = 10 MPa per second, and topological 
plot of the numerical stresses. 
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9^ MPa 



(b) Topological plot (b) Topological plot 


Figure 4.— T ransverse stress concentration, cr-j-j p for an applied Figure 5. — Hydrostatic stress concentration, + a 22 + a 33^* 

stress of =100 MPa at an overall loading rate of for an applied stress of = 100 MPa per second applied at 

1 0 MPa per second after 3600 seconds of creep at 1 00 MPa, an overall loading rate of 6* = 1 0 MPa per second, after 3600 

and topological plot of the numerical stresses. seconds of creep at 1 00 MPa and topological plot of the 

numerical stresses. 
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(a) Hydrostatic stress concentration. 



(b) Topological plot. (b) Topological plot 


Figure 6.— Transverse stress concentration after unloading to Figure 7. — Hydrostatic stress concentration after unloading 

a zero overall stress rate of a® = -10 MPa from the state to a zero overall stress at an overall stress rate of * 

depicted in Fig 4, and topological plot of the numerical -10 MPa from the state depicted in Fig. 5, and topological 

stresses. plot of the numeric stresses. 
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